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We report an extensive and systematic investigation of the multi-point and multi-time correlation functions to 
reveal the spatio-temporal structures of dynamic heterogeneities in glass-forming liquids. Molecular dynamics 
simulations are carried out for the supercooled states of various prototype models of glass-forming liquids 
such as binary Kob-Andersen, Wahnstrom, soft-sphere, and network-forming liquids. While the first three 
models act as fragile liquids exhibiting super- Arrhenius temperature dependence in their relaxation times, 
the last is a strong glass-former exhibiting Arrhenius behavior. First, we quantify the length scale of the 
dynamic heterogeneities utilizing the four-point correlation function. The growth of the dynamic length 
scale with decreasing temperature is characterized by various scaling relations that are analogous to the 
critical phenomena. We also examine how the growth of the length scale depends upon the model employed. 
Second, the four-point correlation function is extended to a three-time correlation function to characterize the 
temporal structures of the dynamic heterogeneities based on our previous studies [K. Kim and S. Saito, Phys. 
Rev. E 79, 060501(R) (2009); J. Chem. Phys. 133, 044511 (2010)]. We provide comprehensive numerical 
results obtained from the three-time correlation function for the above models. From these calculations, we 
examine the time scale of the dynamic heterogeneities and determine the associated lifetime in a consistent 
and systematic way. Our results indicate that the lifetime of the dynamical heterogeneities becomes much 
longer than the a-relaxation time determined from a two-point correlation function in fragile liquids. The 
decoupling between the two time scales is remarkable, particularly in supercooled states, and the time scales 
differ by more than an order of magnitude in a more fragile liquid. In contrast, the lifetime is shorter than 
the a-relaxation time in tetrahedral network-forming strong liquid, even at lower temperatures. 



I. INTRODUCTION 

Various liquids form disordered and amorphous solids 
if temperatures are reduced below their melting points 
while avoiding crystallizations. This transition to a dis- 
ordered solid is known as the glass transitiorPI^. A re- 
markable feature of supercooled states and glasses is the 
drastic increase in the viscosity and the structural relax- 
ation time that accompanies non-exponentially observed 
in various time correlation function^^ltiil. Understanding 
the universal mechanism of the slow dynamics in glass 
transitions is a challenging problem for condensed mat- 
ter physics. 

To tackle this problem, the notion of "spatially hetero- 
geneous dynamics" or "dynamic heterogeneity" in glass- 
forming liquids has attracted much attention in recent 
decades and has been considered central to understand- 
ing the slow dynamics of glasseJ^^^l^. Many theoreti- 
cal, computational, and experimental efforts have been 
devoted to understanding of dynamic heterogeneities in 
glassy systems. 

Experimentally, various nuclear magnetic resonance 
(NMR) and other spectroscopic techniques have revealed 



Electronic mail: kin@ims.ac.jpj 

''^ Electronic mail: .shinji@ims.ac.jp 



"heterogeneous" relaxation of the non-exponential de- 
cays in glassy systems ,^^211211, such systems, the non- 
exponential relaxation can be explained as the superpo- 
sition of indivi dual particle contributions with different 
relaxation rate^^^ED 

A large number of molecular simulations have also pro- 
vided information by allowing for the visualization of mi- 
croscopic details regarding the molecular dynamicd^SMZI, 
These simulations have demonstrated direct evidences of 
the dynamic heterogeneities, i.e., molecular motions ac- 
company correlated domains and, to some extent, exceed 
the molecular scale in a heterogeneous manner in both 
time and space. Experiments to directly visualize these 
dynamic heterogeneities have also been performed in col- 
loidal glasses"*^'^. 

Those results have required characterizing and quan- 
tifying the length and time scales to determine the 
physical role of dynamic heterogeneities in the under- 
lying mechanism of the glassy slow dynamics. Further- 
more, such information regarding spatio-temporal struc- 
tures would be indispensable to an assessment of the 
theoretical scenarios and hypotheses proposed thus far, 
such as the Adam-Gibbs^, random first-order transi- 
tiorP^tm^ dynamic faci litatio rP^, potential en ergy land- 
scap^fil, mode-couplin^SHlS] replicated liquicpSEil^ and 
frustration-limited domain^ approaches. Recent atten- 
tion has focused on determining the physical origin of dy- 
namic heterogeneities by proposing various length scales 
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including those governed by c ooperative rearranging re- 
gions*^^ mosaic lengttP^'SSl, bond-or ientation al order 



(BOOjSHEU^ bond-breakage correlationdSSlsoiill^ ic osahe- 
dral orde'P^^, locally preferred structures (LPSsJ^^'^, 
geometrical frustration22j.SL^ patch correlati ons^^-^ , non- 
local viscoelasticitji^HHl^ Fickian diffusioiPaMIlll, and 
point-to-set (PTS) correlationpHUl. 

Recently, progress has been made in characterizing 
the dynamic length scale via multi-point correlations. 
Extractions of the spatial correlations between particle 
displacements during a typical time in terval ca n lead 
to the four-point correlation functioiPHnHMfini] The- 
oretical treatments have also been provided to analyze 
related multi-poi nt susce ptibilities based on the mode- 
coupling approachlSMinil, Alternative multi-poin t corre- 
lation&i04 i05 loo nUand non-hnear susceptibilitiejnUtllll 
have also been proposed for experimental measurements. 

In addition, a deeper understanding of the time 
scale of heterogeneous dynamics has been provided by 
non-linear responses such as those studies by multi- 
dimensional NMR, hole-burning and photo-bleaching 
techn iques -'^^' ^^ iie-HS,^ and single molecule measure- 
j^gj^^^ii9Hi24| these experiments, the lifetime of the 
dynamic heterogeneity is evaluated from the exchange 
time between the slow- and fast-moving regions, which is 
much longer than the structural relaxation time near the 
glass transition temperature. 

In contrast, much less attention has been paid to theo- 
retical and computational explorations of the character- 
istic time scale of dynamic heterogeneities and its tem- 
perature dependence. Note that a characterization of the 
time scale of dynamic heterogeneities requires an analysis 
of the duration of the heterogeneous motions, which es- 
sentially requires a multi-time extension of the four-point 
correlation function. Some studies have introduced rel- 
evant multi-time correlation functions to investigate the 
time scale of heterogeneitie^^^'^^IIllHIII! However, sev- 
eral time intervals are fixed at a characteristic time scale, 
and thus only limited information is available regarding 
the time scale of the heterogeneous dynamics. 

In our previous 

studie4i2Sli39| 

we have emphasized the 
importance of analyzing a four-point, three-time corre- 
lation function for various time intervals to systemati- 
cally quantify the temporal structures of dynamic het- 
erogeneities and their lifetimes. From the three-time 
correlation function, we have found that the lifetime in- 
creases and becomes much longer than the a-relaxation 
time Tq, of the two-point correlation function for a bi- 
nary soft-sphere supercooled liquid. The observed de- 
coupling between the two time scales at low temper- 
atures is in agreement with the previously mentioned 
experimental results. The exploited strategy is analo- 
gous to the multi-time correlations used in recent multi- 
dime nsional s pectroscopy and other related optical tech- 
nique4^^2H148] These spectroscopic techniques have be- 
come powerful tools for examining the change from het- 
erogeneous dynamics to homogeneous dynamics in vari- 
ous condensed phase systems^^*"^^^. Furthermore, such 



multi-dimensional techniques h ave rece ntly been applied 
to supercooled and glassy stateJ^^^I^^. 

The aim of the present paper is to investigate the 
spatio-temporal structures of dynamic heterogeneities by 
numerically calculating the multi-point and multi-time 
correlation functions. In particular, we extensively in- 
vestigate the manner in which the specifics of the model 
influence the extracted length and time scales. Recently, 
the effects of the model details on the static correlation 
and the dynamics have been critically examined for the 
binary Kob-Andersen model and its Weeks-Chandler- 
Andersen modificatior(i22HilIl While the pair correla- 
tions of these models are almost identical, the structural 
relaxation times differ significantly. These results neces- 
sitate an investigation of the model dependence on the 
length and time scales of the dynamic heterogeneities ex- 
tracted from multi-point and multi-time correlation func- 
tions. Thus, in the present study, we employ a frequently 
used binary mixture of Kob-Andersen, Wahnstrom, soft- 
sphere models, which exhibit a super- Arrhenius tempera- 
ture dependence in their structural relaxation times, and 
a network-forming strong liquid model that exhibits Ar- 
rhenius behavior. We provide comprehensive numerical 
results for the four-point correlation function and the 
three-time correlation function extended from the four- 
point correlation. Using our extensive numerical results, 
we systematically determine the length and time scales 
of dynamic heterogeneities for various glass models. 

The paper is organized as follows. In Sec. [11] we intro- 
duce the simulation details of the glass-forming models 
used in this paper. In Sec. |III[ we present numerical 
calculations of the multi-point and multi-time correla- 
tion functions to characterize the spatio-temporal struc- 
tures of the dynamic heterogeneities 



III A 



we 



In Sec 

briefiy summarize numerical results using conventional 
two-point correlation functions and determine the a- 
relaxation time Tq,. Then, in Sec. |IIIB we evaluate the 



growing length scales of the dynamics heterogeneities in 
the models by analyzing the four-point correlation func- 
tions. Finally, in Sec. |III C[ the lifetimes of the dynamic 
heterogeneities are quantified from the three-time corre- 
lation functions, and their dependence on the model and 
fragility is examined. In Sec. |IV[ we summarize our re- 
sults and provide concluding remarks. 



II. SIMULATIONS OF MODEL GLASSES 

In this study, we carry out extensive molecular dynam- 
ics simulations for three-dimensional binary mixtures in 
the microcanonical ensemble. The system contains A^i 
particles of component 1 and N2 particles of component 
2 under periodic boundary conditions. The total number 
density is fixed at p = N/V with the total number of 
particles N = Ni + N2 and a system volume V. 

The models examined are the well-known proto- 
type models of glass-forming liquids: the binary Kob- 
Andersen Lennard- Jones (KALJ) liquids^^^ the bi- 
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nary Wahnstrom (WAHN) liquidgl^, and binary soft- 
sphere (SS) Uquids^^^ In addition, we also study a 
model of a network-forming (NTW) liqu id th at mimics 
Si02 with short-range spherical potentialJ^^. 



A. KALJ: binary mixture of Kob-Andersen Lennard-Jones 
particles 

The binary Lennard-Jones mixture is the most fre- 
que ntly uti lized model for the study of glass transi- 
tionJi^lIlM]^ The pair potentials are given by 



potentials 



(1) 



in which a,/3 S {1,2} are the indexes of the particle 
species. The energy and size ratios are ei2/eii = 1-5, 
£22/611 = 0.5 and cri2/o'ii = 0.8, (T22/fii = 0.88, re- 
spectively. The masses of the two particle species are 
equal, mi = m2 = 1. The interaction is truncated 
at r = 2.5crQ^. The reduced units an, en/ks, and 
\/miali/eii are used in the model for length, tempera- 
ture, and time, respectively. The time step is At — 0.001 
in the reduced time units. The total number density 
is fixed at p = 1.2 with Ni = 800 and 7V2 = 200. The 
temperatures investigated are T = 0.7, 0.65, 0.6, 0.55, 0.5, 
and 0.47. 



(3) 



with the cubic smoothing function Va/sir) — B{a~r)^ +C 
for distances r > Tc = \/3- The values of a, B, and C are 
determined by the continuity conditions up to the second 
derivative of Vap (r). The size, mass, and energy ratios are 
the same as those of the WAHN model: Oxjai — 1/1.2, 
mi/m2 — 1/2, and ei/e2 = 1, respectively. Thus, this 
model can be regarded as a purely repulsive interact- 
ing system of the WAHN model. Simulation results will 
be described in terms of the reduced units cti, ti/ks, 
and \pm\a\J(~\ for length, temperature, and time, re- 
spectively. 

The thermodynamic state of this model is usually char- 
acterized by the following non-dimensional parameter: 



P 



ei 



1/4 



7 3 



(4) 



in which Zq represents the effective particle size defined 
by 4^0^ = (2cri)3 + 2(cri ^ 02)^ ^ {'^oif . In the sim- 
ulation, the total number density is given as p = 
with Ni = N2 = 500. The investigated states are 
r = 1.30, 1.35, 1.38, 1.42, 1.45, and 1.47. The correspond- 
ing temperatures are T = 0.350, 0.301, 0.276, 0.246, 0.226, 
and 0.214 with a time step of At 0.005. 



B. WAHN: binary mixture of Wahnstrom Lennard-Jones 
particles 

The KALJ model is a non-additive mixture and thus 
disobeys the so-called Lorentz-Berthelot combining rules 
due to the strong attraction between components 1 and 
2. Alternatively, the prototype model of the additive and 
equimolar bin ary Len nard-Jones mixture is introduced 
by WahnstronP^mm^ Xhe interaction potentials are the 
same as in Eq. ([T]), in which the size, mass, and energy 
ratios are given as uxjui = 1/1-2, TOi/m2 = 1/2, and 
£1/62 = 1, respectively. The Lorentz-Berthelot rules. 



(2) 



are obeyed in this model. Simulation results will be de- 
scribed in terms of the reduced units cti, ei/fcs, and 
mial / ei for length, temperature, and time, respec- 
tively. The system consists of iVi — 500 and N2 = 500 
particles with a fixed number density p = 0.75. A time 
step of At = 0.001 is used. The temperatures investi- 
gated are T = 0.8, 0.75, 0.7, 0.65, 0.6, and 0.58. 



C. SS: binary mixture of soft-sphere particles 

We also stu dy an equimolar binary mixture of soft- 
sphere particlesISSIlSZl, Particles interact via the soft-core 



D. NTW: tetrahedral network-forming liquids 

In addition, we study a model of network- forming liq- 
uids interacting via spherical short-ranged potentials'^. 
This model is simple model and imitates Si02 glasses, 
in which tetrahedral networks strongly dominate the dy- 
namics with an Arrhenius behavior for the structural re- 
laxation time, even near the glass transition temperature. 
The interaction potentials are given as 



r 



12 



(1 " ^o.,) ("^ 



(5) 



in which Sap is the Kronecker delta. The interaction is 
truncated at r = 2.5crQ^. The size and energy units are 
determined as follows: 



O'12/o'll 

ei2/eii 



0.49 
24 



O'22/fll 
£22/ £11 



0.85 
1. 



(6) 
(7) 



The mass ratio is determined as m2/mi — 0.57 from the 
same ratio of O and Si. The units of length, time, and 
temperature are given as cth, en/ks, and \/m]a\^[[eii, 
respectively. These parameters are adjusted to reproduce 
the radial distribution functions of the Si02 amorphous 
states. Tetrahedral networks are found to be formed due 
to the highly asymmetric size ratio and the strong at- 
traction between the different components'^*. The den- 
sity of the investigated system is p = 1.655, with parti- 
cle numbers Ni = 1,000 and N2 = 2,000. This value 
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FIG. 1. Self-part of the intermediate scattering function 
Fsik, t) of the component 1 particles for various glass-forming 
liquid models: (a) KALJ, (b) WAHN, (c) SS, and (d) NTW. 



corresponds to the density p = 2.37gA~'^ of the so- 
called van Be est-Kra mer-van Santen (BKS) model for 
the silica glasJ^^^HIZll, xhe simulations are carried out at 
T = 0.42, 0.4, 0.38, 0.36, 0.34 and 0.32 with a time step of 
M = 0.0005. 



III. RESULTS AND DISCUSSION 

A. Intermediate scattering function and the a-relaxation 
time 

First, we study the conventional two-point density 
correlation function and determine the structural a- 
relaxation time Tq,. The self-part of the intermediate 
scattering function of the component 1 particles 



(8) 



i=i 



is calculated for various glass-forming models. Here, 
Arj(0, t) = rj{t) — rj{0) is the jth particle displacement 
vector at times and t, and k is the wave vector. 

The behavior of Fs(k, t) at various temperatures is 
demonstrated in Fig. IT] The wave vector k — \k\ is 
chosen so that the static structure factor of component 
1, S'ii(fc), marks the main peak as (a) k — 7.25, (b) 
6.65, (c) 6.55, and (d) 8.0 for the KALJ, WAHN, SS, and 
NTW models, respectively. From these calculations, we 
determine the a-relaxation time Tq as Fs{k,Ta) = e~^. 
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FIG. 2. (a) Qf-relaxation time Tq as a function of the inverse 
temperature To/T with the onset temperature Tq. (b) q- 
relaxation time Ta as a function of the temperature difference 
(T — Tc)/Tc from the mode-coupling divergence temperature 
Tc- Each dashed line represents the power-law fits with ex- 
ponent A = 2.4, 1.8, 2.2, and 2.8 for the KALJ, WAHN, SS, 
and NTW models, respectively. 



In Fig. [2][a), the temperature dependence of Ta is plot- 
ted as a function of the inverse temperature Tq/T. Here 
Tq denot es t he onset temperature introduced in the ear- 
lier worlPfil. We set To as Tq = 0.8, 0.7, 0.3, and 0.5 
for the KALJ, WAHN, SS, and NTW models, respec- 
tively. Below this onset temperature, Fs{k,t) begins to 
develop a two-step relaxation in each model. Further- 
more, the temperature dependence of Tq, exhibits super- 
Arrhenius behavior in the KALJ, WAHN, and SS models. 
This behavior is typical for the fragile glass-forming liq- 
uids. In contrast, in the NTW model the tetrahedral net- 
works begin to develop strongly below Tq. Correspond- 
ingly, the temperature dependence of the a-relaxation 
time obeys the Arrhenius law^^^. In addition, we sum- 
marize the power-law behavior as Tq. cx IT — Tc\~^, as 
predicted by the mode-coupling theorjlSU£2l The re- 
sults of the power-law fittings for the four model liq- 
uids are demonstrated in Fig.[2]^b). We obtain the values 
of the exponent A and the mode-coupling temperature 
Tc as (A,Tc) « (2.4, 0.435), (1.8, 0.56), (2.2, 0.198), and 
(2.8,0.31) for KALJ, WAHN, SS, and NTW liquids, re- 
spectively. Similar results have been reported for the 
KALJ, SS, and NTW modelsPal. We also note that the 
power-law behavior of the NTW model is reliable over the 
limited temperature ran ge, a s investigated in the simula- 
tions for the BKS modeP^Sl. 



B. Four-point correlations and the growing length scale of 
dynamic heterogeneities 

To characterize the growth of the dynamic hetero- 
geneities in supercooled states, the four-point correlation 
function is introduced to measure the correlation of the 
particle mobility field at a given time interval. There are 
several choices for the mobilitj^field such as the the par- 
ticle displacement amplitude^'^', the overlap function'^, 
and the intermediate scattering functiorPH. Physically, 



the choice of the mobiHty field does not alter the fun- 
damental meaning of the four-point correlation function. 
We here calculate the four-point dynamical susceptibility 
X4{k,t), which is defined as the variance of the self-part 
of the intermediate scattering function Fs{k,t): 



(a) KALJ 



(b) WAHN 



We utilize Fs{k,t) expressed as 



1 

Fs{k,t) = — ^cos[fc. Ar,(0,t)], 



(9) 



(10) 



i=i 



with Fs{k,t) = (Fs(fc,t)). The total value of Xiik,t) is 
approximately proportional to the extension of the spa- 
tial correlations in dynamics with fc at a given time inter- 
val t because Xiik, t) investigates the increasing deviation 
of the two-point correlation function Fs{k,t) from the 
mean behavior, As is well documented in various stud- 
ies^"^"^'^'' and as demonstrated in Fig.jsj X'i{k, t) typically 
has its maximum value at the time scale of Tq, which 
increases as the temperature decreases. The growth of 
X4(fc,i) for the strong NTW hquid is suppressed and is 
smaller than those of other fragile liquids. This behav- 
ior implies that the dynamic heterogeneity is less pro- 
nounced and plays a mino r role in th e strong liquid, as 
revealed in previous 

studiefpSEHEZIl. Therefore, the dy- 
namics of the strong liquid can be interpreted as mostly 
occurring through the rearrangement of strongly connect- 
ing tetrahedral networks. Interestingly, a similar sup- 
pression of X4:{k, t) in strong liquids, for l onger times, has 
bee n obse rved in polydispersed systemd^^^E^, c olloidal 
ggjJmHH!^ and confined systems in random medicP^^^l^. 

The study of the spatial correlation of the four-point 
correlation function is also essential to extract the grow- 
ing length scales ^ of the dynamic heterogeneities. In this 
study, instead of the self-intermediate scattering func- 
tion, we utilize the frequently used four-point correlation 
function using the overlap function, which is defined as 
follows: 



S4{q,t) = ^{Q{q,t)Q{-q,t)), 



Qiq,t) 



^ VFj(a,t) exp[ 



iq ■ r 



(11) 



(12) 



with the wave vector q 
e(a - \r,it) - r, 

the Heaviside step's function 0(a:'3SH2l 



\rj{t) — rj(O)l) is the overlap 



Here, Wj{a,t) — 
ing function with 
The function 



Wj[a,t) selects a particle that moves further than the 
distance a during the time interval t. The value a = 
0.3 is typically chosen. As studied in the previous 
studjS^, the relaxation profile of the overlap function 

Q{t) = (lM)(E^i W^i(a,i)) with this probe length 
scale a = 0.3 approximately corresponds to that of 
the self-intermediate scattering function Fs{k, t) with the 
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FIG. 3. Four-point dynamical susceptibility Xi{k,t) for the 
(a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. 



wave number marking the main peak of the static struc- 
ture factor. We also note that any significant differ- 
ences in Si{q,t) are not observed if we choose Fs{k,t) 
as the mobility field. Although recent reports have de- 
scri bed the res ults of Si{q,t) using large-scale simula- 
tion d^^^ l l ^Utl^, we revisit the identification of the corre- 
lation length ^ from the four-point correlation function, 
Eq. (11). Accordingly, we use much larger systems with 
N = 100, 000 for the KALJ, WAHN and SS liquids and 
N = 90,000 for the NTW liquid to calculate Si{q,t). 



The other parameters described in Sec.|lljare unchanged 
throughout the simulations. Correspondingly, the linear 
dimension of the system, in termed of the unit length, 
L = V^/^, is given as i = 43.68, 51.09, 51.27, and 37.89 
for the KALJ, WAHN, SS, and NTW models, respec- 
tively. 

Figure |4] shows the wave- number dependence of S/i{q, t) 
on the time scale r^. As indicated in Fig. |4j Si{q,t) at 
the a-relaxation time grows in small wave-numbers of 
q, particularly at low temperatures. This observation 
indicates that the mobile (or immobile) particles become 
highly correlated in space when the system undergoes 
supercooling. The behavior of S^^qyt) at small wave- 
numbers can be described by the Ornstein-Zernike (OZ) 
form: 



Si{q,t) = 



xit) 



(13) 



in which ^(t) is regarded as the length scale of the dy- 
namic heterogeneity and x(t) is the dynamic suscepti- 
bility at g — 0. The OZ form with a = 2 has been 
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FIG. 4. Four-point static structure factor S4,{q, To) at various 
temperatures as a function of the wave-number q for tlie (a) 
KALJ, (b) WAHN, (c) SS, and (d) NTW models. 



used in previous studiej ^ * '^^ * ^^ [ However, Fig. [4] shows 
that the exponent a depends on the details of the sim- 
ulation model. Figure [5] displays the scaled function 
SiiqiTa) Ixija) as a function of qiija)- The results are 
in good agreement with a ~ 2.4, 2.0, 2.4, and 1.5 for 
the KALJ, WAHN, SS, and NTW models, respectively. 
The apparent difference in the exponent a between the 
fragile (KALJ, WAHN, and SS) and strong (NTW) glass- 
forming liquids may be related to the change in geometri- 
cal characteristics of the heterogeneous mo tion s if we em- 
ploy an analogy to the critical phenomenaP^. A similar 
difference in the exponent has been found in kinetically 
constrained models (KCMs), in which a snapshot of the 
dynamic heterogeneity in strong KCM model exhibits a 
smoother cluster structure than in fragile KCM modePi'. 

The values of the length scales ^ and the susceptibility 
X at the a-relaxation time are determined from the fitting 
to Eq. (13). Figure [6](a) shows the temperature depen- 
dence of the correlation length ^(tq ) . We find that the in- 
creasing length scale ^(tq) with decreasing temperature 
can be described by the mode-coupling-like power-law 
^{to) ~ \T—Tc\~'' at the investigated temperatures. The 
exponent is w 0.5 for the fragile KALJ, WAHN, and SS 
models, whereas v decreases to v ss 0.25 for the strong 
NTW model. However, the length ^ at lowest tempera- 
ture T — 0.32 is apparently deviated from the power-law 
behavior of the NTW model, which is a strong liquid ex- 
hibiting the Arrhenius behavior. This limited power-law 
behavior is also observed in the temperature dependence 
of the a-relaxation time (see Fig. [2|^b)). In addition. 



FIG. 5. Scaled four-point static structure factor 

54(9, Ta)/x(''"a) as a function of q(,{Ta) at various tempera- 
tures for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW 
models. The dashed line represents the Ornstein-Zernike 
form 1/(1 -I- (q^ir^))") with (a) a = 2.4, (b) 2.0, (c) 2.4, 
and (d) 1.5. 



we obtained the scaling relationship Tq ^ ^{tcY found 
in Fig. |6f^b). In analogy to the dynamical critical phe- 
nomena, the exponent z is approximately determined by 
the relationship z — A/v, i.e., the exponent z becomes 
z w 4.4, 4.8, 3.6, and 11.0 for the KALJ, WAHN, SS, and 
NTW models, respectively. 

We also examine the relationship between the suscepti- 
bility x(''q) and the correlation length ^(tq,) via the scal- 
ing xiTa) ~ ^i'''a)'^~^7 as demonstrated in Figure |6][d) 
These results indicate that the scaling exponent of each 
model is correlated to the value of the exponent a utilized 
in the OZ function, as shown in Eq. (13 1. In addition. 



the dynamic susceptibility xi^a) is approximated by the 
scaling relation x('^q) ~ |T — Td"'^ with 7 = (2 — r])^, as 
shown in Fig. [6jc) . Similar dynamic criticality has been 
investigated in Refs. BTIand llOOl 

Here, we remark that a recent theoretical treatment 
has been presented based on mode-coupling approach, 
referred to as the inhomogeneous mode-coupling theory 
(IMCT)|i2Ili2a). The IMCT predicts the scaling exponents 
to be i/ 1/4 and 2 — 77 = 4 at the a- relaxation time 
scale. To examine the validity of these predictions, in- 
tensive large-sc ale molecular simulations have been car- 
ried out!i^i"i^. Our exponents, z « 4.4 and 2 — 77 « 2.4, 
of the KALJ model are similar to the values obtained in 
Refs. [T04l and [T82l In these findings, the molecular sim- 
ulations disagree with the theoretical predictions. How- 
ever, the IMCT does not investigate a four-point correla- 
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FIG. 6. (a) Correlation length ^(tq) as a function of 
{T — Tc)/Tc. The dashed line represents the power-law behav- 
ior, ~ |T - Tel"" with u = 0.5 (red) and f = 0.25 (or- 
ange), respectively, (b) The relationship between the correla- 
tion length ^{tci) and the a-relaxation time Tc- The dashed 
line represents the power-law behavior, Tc ~ ^{tc)^ with z = 
4.4 (red) and z = 11.0 (orange), respectively, (c) The dynamic 
susceptibility x{'''a) as a function of {T — Tc)/Tc. The dashed 
line represents the power-law behavior, x(''"a) ~ |T — Tel"'' 
with 1/ = 1.2 (red) and f = 0.37 (orange), respectively, (d) 
The relationship between the correlation length ^{tc) and the 
dynamic susceptibility x{'''a)- The dashed line represents the 
power-law behavior, x('''q) ~ £,{Ta)^^^ with 2 — rj — 2.4 (red) 
and 2 ~ rj — 1.5 (orange), respectively. 



tion function such as Eq. (11) but treats the three-point 



correlation defined as the response function of the two- 
point correlation function due to the inhomogeneous ex- 
ternal field. Further analysis is needed to understand the 
relationship between the four-point correlations utiHzed 
here and the inhomogeneous three-point susceptibilities. 



C. Three-time correlations and lifetimes of dynamic 
heterogeneities 

In this subsection, we explore and accentuate the char- 
acteristic time scale of the dynamic heterogeneity. To 
this end, we search for the multi-time extension of the 
four-point correlation function Xiikjt) by introducing 



/ 1 



3 = i 



(14) 
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FIG. 7. Two-dimensional correlation maps of the three-time 
correlation functions AF4(fc, ti, ^a) for the (a) KALJ, (b) 
WAHN, (c) SS, and (d) NTW models. The state is chosen 
at the lowest temperature for each model as (a) T = 0.47, 
(b) T = 0.58, (c) r = 1.47, and (d) T = 0.32. The waiting 
times t2 normalized by Tq are increased from left to right in 
each model. The vertical (horizontal) dashed line denotes the 
a-relaxation time as ti = Tq, (ts = Tq). Note that the waiting 
time is different in each figure. 



in which 



cos[k ■ Arj(0,t)] - Fs{k,t) 



(15) 



provides the individual fluctuations in the two-point cor- 
relation function at times and t. This three-time cor- 
relation function examines the correlations at four differ- 
ent times, 0, ri, T2, and T3. In practice, Ai^4(fc, ii, ^2, ^3) 
reveals the correlation of the two-point correlation func- 
tion Fg{k,t) between the two time intervals, ti — ti and 
^3 = — T2. Furthermore, the progressive changes in 
the waiting time t2 = T2 — ti of the three-time corre- 
lation function A_F4(fc, ti, ^3) allow for an investiga- 
tion of how the correlated motions change with time 
t2. In fact, the first time-interval portion SFj{k,0,Ti) 
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FIG. 8. The diagonal portions of the three-time correlation 
functions AF4{k,t,t2,t) for the (a) KALJ, (b) WAHN, (c) 
SS, and (d) NTW models. The state is chosen at the lowest 
temperature for each model as (a) T — 0.47, (b) T = 0.58, 
(c) r = 1.47, and (d) T = 0.32. The waiting times t2 are 
normalized by Tc at each temperature. 



enables a distinction between the subensemble of slow- 
and fast- moving particles in the dynamic heterogeneities. 
In addition, the total function 6Fj{k,0,Ti)6Fj{k,T2,T3) 
can reveal how long the correlations in the dyna mics of 
the subensembles remain over the waiting time ^ ji38 | i39 | 
Therefore, this multi-time correlation function can pro- 
vide the temporal structures and the associated charac- 
teristic time scale of the dynamic heterogeneity. In other 
words, the time scale extracted from the three-time cor- 
relation function can be regarded as the lifetime of the 
dynamic heterogeneity, which should be associated with 
the length scale ^ determined from the four-point corre- 
lation function. 

Figure [T] shows the time evolutions of the three-time 
correlation functions AF4 at the lowest temperature for 
each glass-forming liquid. Note that the wave-number 
k is chosen as the same value used in the calculation of 
Fs{k,t) in Fig. [1] We also provide the diagonal portions 
of the time evolutions along t = ti = for various t2 
in Fig. [sj In previous studie^^^lIIMl^ have examined 
the three-time correlation functions for the SS model. 
This work confirms that the basic features of AF4 are 
similar and that the time evolutions at ^2 occur similarly 
among all of the glass models studied: (i) As the tempera- 
ture decreases, the intensity of the three-time correlation 
AFi{k,ti,t2,t3) gradually increases (see the result s fo r 
other temperatures in the Supplementary MateriaPSll) . 
This result is due to the correlations between particles 



that move slower (faster) during the first time interval ti 
and remain slower (faster) during the second time inter- 
val t^. (ii) The correlations of AF4(A;, <i, ^2, ^3) aX t2 = 
between the first time interval ti and the subsequent time 
interval are noticeable over wide time scales. This re- 
sult implies that the particle motions are coupled not 
only at the (diagonal) a- and a-relaxation time scales, 
but also at the (off-diagonal) a- and /3-relaxation time 
scales, (iii) With increasing the waiting time t2, AF4 
gradually decays to zero, indicating that the dynamics 
change from heterogeneous to homogeneous because of 
the memory loss in the correlated motions between the 
two time intervals ti and ts for the given waiting time 
t2- However, it should be noted that the lineshape of the 
NTW model, which exhibits the strong (Arrhenius) glass 
behavior, differs from those of the KALJ, WAHN, and SS 
models, which exhibit fragile glass (super- Arrhenius) be- 
havior, particularly at t2 = 0. We observe the strong 
correlations of AF4 at small values of ti and in the 
NTW model, which approximately corresponds to the 
time scale of the early /3-relaxation at which Fs{k,t) 
exhibits damped oscillations on a plateau as shown in 
Fig. flVd). T his well-known finite-size effect in silica 
glasseP^^Em would be smaller if simulations were car- 
ried out for larger systems. Furthermore, the relaxation 
time scale of AF4(fc, ii, ^2, ^3) depends on the model. As 
demonstrated in Fig. ^d), AF4 of the NTW model de- 
cays rapidly. This time scale is clearly smaller than the a- 
relaxation time determined by the two-point correlation 
function. In contrast, the relaxation of A_F'4 occurs on a 
time scale comparable to (or exceeding) Tq in the fragile 
KALJ, WAHN, and SS models, as depicted in Fig. [Wa)- 
(c). 

To quantitatively distinguish between the models we 
quantify the average lifetime of the dynamic hetero- 
geneities using the waiting-time t2 dependence of A_F'4. 
To this end, we define the volume of the heterogeneities 
by integrating over the two time intervals ti and tgiSHIlS!^ 

/>oo poo 

Ahctcro(i2)= / dt3 dHAFiikMMM)- (16) 
Jo Jo 

This integration resembles the underlying strategy of 
non-linear responses such as NMR, ho le-burning, and 
photo-bleaching technique In simulations, 
similar procedures hav e been utilized to analyze relevant 
multi-time correlationji^llilll. Figure [9] illustrates the ^2 
dependence of AhctGro(i2) normalized by Ahetero(O) for 
the KALJ, WAHN, SS, and NTW models. The waiting 
time t2 is normalized by at each temperature. Fig- 
ure |9][d), clearly shows that Ahctcro of the NTW model 
decays rapidly and that the relaxation time is smaller 
than Ta at any temperature. In contrast. Fig. |9ja)-(c) 
demonstrates that in the fragile liquid models, the relax- 
ation of Ahctcro is slower than Tq, with decreasing tem- 
perature. Remarkably, the characteristic time scale of 
Ahetero for the WAHN model at the lowest temperature 
exceeds t„ by more than an order of magnitude. 

The dependence of the normalized volume 
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FIG. 9. Waiting time t2 dependence of the integrated three- 
time correlation function Ahetero(fc, fa)/Ahetero(fc, 0) for the 
(a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The 
waiting times are normalized by Tc for each temperature. 
The solid curve is determined by a fitting with the stretched- 
exponential form exp[— (f2/Thotoro)'^] with c « 0.6, 0.5, 0.5, 
and 0.3 for the KALJ, WAHN, SS, and NTW models, respec- 
tively. 
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FIG. 10. (a) Average lifetime of DH rhotoro normalized by the 
a-relaxation as a function of the inverse temperature Tq/T 
with the onset temperature Tq. The dashed line represents 
Thetero = Tq as a vlewiug guide, (b) The relationship between 
two time scales, Thotcro and Ta. The dashed line is the power- 
law behavior Thetero ~ tJ^ with a slope of « 1.9, 1.5, 1.25, 
and 0.9 from top to bottom. 



the dynamic heterogeneity as Thctoro at various temper- 
atures for each model. We plot Thotcro as a function of 
the inverse temperature Tq/T in Fig. [lO{^a). As shown 
in Fig. [To{^b), we obtain a relationship between the two 
time scales Thotcro and Tq, that follows the power-law-like 
behavior, Thotcro ^ Tq'' with Q « 1.25, 1.9, 1.5, and 0.9 
for the KALJ, WAHN, SS, and NTW models, respec- 
tively. From the analysis, we find that Thotcro of the 
network-forming strong glass (NTW) is not greater than 
Ta and tends to decrease with decreasing temperature T. 
That relationship is responsible for the minor role of the 
dynamic heterogeneities in strong liquids, as previously 
discussed when we examined the four-point correlations 
and the associated length scale ^ in Sec. |IIIB[ In 
contrast, the ratio Thotcro /''"q in the fragile liquid models 
exhibits the opposite temperature dependence, i.e., 
the lifetime Thotcro exceeds the a-relaxation time Tq, 
with decreasing temperature. However, the increase 
in Thotcro upon supercooling is considerably different 
among the fragile KALJ, WAHN, and SS models. The 
ratio Thotcro/Ta markedly exceeds the unity at lower 
temperatures in the WAHN and SS models, whereas 
Thetero m the KALJ model remains on the time scale of 
Ta, even at the lowest temperature. We note that this 
feature of the KALJ model was also found in a previous 
study using a similar multi-time corr elat ion function, in 
which its lifetime is comparable to TqP^. 

A discussion of the significant observed dependence of 
the lifetime Thotcro on the model details is meaningful, 
even for the fragile models. Recently, the fragility in- 
dexes of the simulated models such as the KALJ and 
WAHN models have been critically investigated from the 
perspective of the many-body static correlations hidden 
in the two-point correlations such as the usual radial dis- 
tribution function. It has been found that the slow- and 
long-lived correlated domains correspond to the locally 
preferred structur es (LP Ss) , that are characterized by the 
Voronoi polyhedrgP^Hl. These studies have also revealed 
that the non-additive KALJ mixture is less fragile than 
the additive WAHN mixture. This difference in fragility 
can be explained in terms of the spatial extent of the 
LPS domains. In fact, the growth of the LPS domains is 
significant in the WAHN model that develops icosahedral 
order upon supercooling. In contrast, the LPS domains 
in the KALJ model formed by a bicapped prismatic order 
are found to be smaller than that in the WAHN model. 
Given these findings, one can conclude that the over- 
all temperature tendency of the lifetime Thotcro shown in 
Fig.[TO]is correlated and intensely sensitive to the fragility 
of the model, i.e., more fragile liquids tend to exhibit 
longer dynamic heterogeneity lifetimes Thctoro- 



Ahotcro(i2)/Ahctero(0) ou the waiting time is can 
be fitted by the stretched-exponential function 
exp[—(t2 /Thotcro)'^], as demonstrated in Fig. [9j The 
exponent c is approximately c ~ 0.6, 0.5, 0.5, and 0.3 for 
the KALJ, WAHN, SS, and NTW models, respectively 
From this analysis, we determine the average lifetime of 



IV. SUMMARY 

We have examined the four-point correlation function 
and its three-time correlation extension to systematically 
characterize the length and time scales of dynamic het- 
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erogeneities in prototype fragile (KALJ, WAHN, and SS) 
and strong (NTW) glass models. Analyses such as the ex- 
tensive investigation performed herein to determine not 
only the length scale but also the time scale of various 
glass- forming models have not, to the best of our knowl- 
edge, been previously reported. 

First, we quantified the growing length scale of the dy- 
namic heterogeneities upon supercooling as determined 
by the wave-number dependence of the four-point corre- 
lation function using the overlap function. The scaling 
relationships of the extracted length scale which are 
analogous to the dynamical scaling obtained for critical 
phenomena, were consistently explored for the employed 
glass models. We observed that the length scale increases 
with decreasing temperature depending on the fragility 
of glass. In particular, the increase in the dynamic length 
scale of the strong glass is suppressed compared with 
those of the fragile glass-forming liquids, indicating that 
the dynamic heterogeneity is less pronounced and plays 
only a minor role in the strong liquid. We also com- 
mented on the comparisons of our numerical results with 
the theoretical predictions of the IMCT. 

Second, we investigated the time scale of the dynamic 
heterogeneities from determining how long the heteroge- 
neous dynamics survive. Comprehensive numerical re- 
sults of the three-time correlation function were demon- 
strated via two-dimensional correlation maps with an 
analogy to the multi-dimensional spectroscopic meth- 
ods, as outlined in the introduction. From the progres- 
sive changes in the second time interval of the three- 
time correlation function, we quantified the character- 
istic time scale of the dynamic heterogeneities and the 
associated lifetime Thetero- The lifetime Thotcro exceeds 
the a-relaxation time Tq, particularly for highly super- 
cooled states in fragile glass-forming liquids. In contrast, 
Thetero IS Smaller than Tq, even at lower temperatures in 
the strong liquid, indicating that the dynamic hetero- 
geneities play a minor role. Furthermore, we observed 
that the temperature dependence of Thetero depends sig- 
nificantly on the fragility, i.e., more fragile liquids exhibit 
long-lived dynamic heterogeneities with a time scale that 
exceeds the a-relaxation time. The two time scales dif- 
fer by more than an order of magnitude in the WAHN 
model. 

Finally, we remark that it is of important to investi- 
gate the relationship between the length and time scales 
of dynamic heterogeneities and their model dependence. 
In fact, we find that the length scale ^ of fragile liquids 
increases with decreasing temperature in a similar man- 
ner, as seen in Fig.[6]ja). In contrast, the time scale Thotcro 
is more sensitive to the fragility and becomes noticeably 
longer than t^ as the fragility index increases. To clar- 
ify it, further work that utilizes not only multi-point and 
multi-time correlations but also other measurements in- 
cluding configurational entropy, LPS, PTS, and BOO is 
required. 
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